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Abstract.  The  multi-domain  hybrid  Spectral- WENO  method  (Hybrid)  is  introduced  for  the 
numerical  solution  of  two  dimensional  nonlinear  hyperbolic  systems  in  a  Cartesian  physical 
domain  which  is  partitioned  into  a  grid  of  rectangular  subdomains.  The  main  idea  of  the  Hy¬ 
brid  scheme  is  to  conjugate  the  spectral  and  WENO  methods  for  solving  problems  with  shock 
or  high  gradients  such  that  the  scheme  adapts  its  solver  spatially  and  temporally  depending 
on  the  smoothness  of  the  solution  in  a  given  subdomain.  Built  as  a  multi-domain  method, 
an  adaptive  algorithm  is  used  to  keep  the  solutions  parts  exhibiting  high  gradients  and  dis¬ 
continuities  always  inside  WENO  subdomains  while  the  smooth  parts  of  the  solution  are  kept 
inside  a  spectral  one,  avoiding  oscillations  related  to  the  well-known  Gibbs  phenomenon  and 
increasing  the  numerical  efhciency  of  the  overall  scheme.  A  higher  order  version  of  the  multi¬ 
resolution  analysis  proposed  by  Harten  is  used  to  determine  the  smoothness  of  the  solution 
in  each  subdomain.  We  also  discuss  interface  conditions  for  the  two  dimensional  problem  and 
the  switching  procedure  between  WENO  and  spectral  subdomains.  The  Hybrid  method  is  ap¬ 
plied  to  the  two-dimensional  Shock- Vortex  Interaction  and  the  Richtmyer-Meshkov  Instability 
(RMI)  problems. 
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1  Introduction 

In  this  article  we  extend  the  one-dimensional  Multi-Domain  Hybrid  Spectral- WENO  Method 
(Hybrid)  for  Hyperbolic  Conservation  Laws  [6]  to  two-dimensions  in  space  and  apply  it  to  the 
classical  Shock- Vortex  interaction  and  Richtmyer-Meshkov  instabilities  problem.  The  general 
idea  of  the  Hybrid  method  is  to  use  a  multi-domain  framework  in  order  to  apply  convenient 
spatial  discretizations  to  the  smooth  and  rough  parts  of  the  numerical  solution.  Shocks  and 
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high  gradients  are  kept  at  WENO  subdomains,  while  complex,  but  still  smooth,  details  of  the 
solution  are  treated  within  spectral  subdomains.  Numerical  efficiency  is  increased  with  respect 
to  the  classical  spectral  and  WENO  methods:  Postprocessing  techniques  of  the  spectral  method 
approach  of  shocks  [11,  16]  are  avoided,  since  no  Gibbs  phenomenon  will  occur,  and  the  expensive 
characteristic  decompositions  and  projections  of  the  WENO  method  are  skipped  at  the  smooth 
parts  of  the  solution  [2,  3,  5,  9,  10,  15]. 

The  main  issues  in  the  construction  of  the  Hybrid  method  are  the  smoothness  measurement 
of  the  solution  and  the  subdomains  types  switching  algorithm.  In  this  work  we  employ  the 
high  order  multi-resolution  algorithm  by  Ami  Harten  [12]  to  build  a  local  classification  of  the 
solution  into  smooth  and  rough.  Originally  built  to  decrease  the  work  of  the  fluxes  computations 
of  Conservation  Laws,  Harten’s  Algorithm  proposes  to  use  information  from  coarser  grids  when 
the  solution  is  locally  over-represented.  We  instead  use  the  multi-resolution  information  to  apply 
distinct  numerical  methodologies  to  the  different  structures  of  the  solution.  The  main  goal  is 
to  conjugate  the  higher  efficiency  of  the  spectral  method  with  the  shock-capturing  capability  of 
the  WENO  method.  The  multi-resolution  analysis  is  used  to  trigger  the  switching  algorithm  to 
change  the  subdomains  spatial  discretizations  if  shocks  start  to  develop  at  a  spectral  subdomain, 
or  if  the  solution  becomes  smooth  at  a  WENO  one.  Moving  discontinuities  are  similarly  treated 
by  changing  to  (or  maintaining  as)  WENO  the  subdomains  on  their  paths  and  switching  to  (or 
maintaining  as)  spectral  the  subdomains  that  were  left  behind.  These  changes  are  performed 
via  Lagrangian  and  spectral  interpolations  of  the  local  solutions  to  the  new  discretizations  grids. 
Interpolation  is  also  used  to  patch  the  solutions  at  the  interfaces.  While  a  simple  average  is 
sufficient  for  the  interfaces  where  the  solution  is  smooth,  using  the  same  grid  spacing  at  adjacent 
WENO  subdomains  is  necessary  for  a  conservative  transmission  of  shocks  [24].  Even  though  we 
do  not  have  a  theoretical  proof  of  the  conservation  of  the  Hybrid  scheme,  we  argue  that  with 
the  conjugation  of  two  conservative  schemes  with  a  conservative  WENO  interface  and  the  high 
order  accuracy  of  the  conservative  spectral  scheme,  the  conservation  error  should  be  spectrally 
small.  We  have  numerically  demonstrated  this  fact  in  [6]  through  a  long  time  integration  of  the 
inviscid  Burgers  equations  with  correct  shock  speed  and  achieved  excellent  agreement  with  the 
analytical  solution  of  the  standard  Riemann  shock-tube  problems,  such  as  the  Lax  and  the  Sod 
problem  of  the  Euler  Equations. 

The  paper  is  organized  as  follows:  Section  2  provides  quick  reviews  on  spectral  and  WENO 
methods.  The  Multi-Resolution  analysis  is  discussed  in  details  at  Section  3  and  the  Hybrid 
Method  is  introduced  at  Section  4.  The  Switching  Algorithm  is  presented  in  Section  4.3  and 
numerical  experiments  with  two  dimensional  compressible  flows  are  finally  presented  at  Section 
5.  Concluding  remarks  are  given  in  Section  6. 

2  Spectral  and  WENO  Methods 

2.1  The  Spectral  Method 

In  spectral  collocation  methods,  the  function  u{x)  is  interpolated  by  a  global  Lagrangian  inter¬ 
polation  polynomial  of  degree  N,lj{x),  at  a  given  set  of  collocation  points  {xj,j  =  0,...,N} 
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as 


N 

Imu{x)  =  E  u{xj)lj{x)  , 

j=0 


where  1^  is  the  Interpolating  operator  and  lj{xi)  =  5ij. 

The  Lagrangian  interpolation  polynomial  lj{x)  can  be  constructed  as 


(2.1) 


lj{x) 


q{x) 

{x  —  Xj)q'{xj)  ’ 


and  the  derivative  of  I^u  becomes 


N 

q(.x)  =  JJ(x  -  Xj) 
j=0 


dlffu{x) 

dx 


N 

'^u{Xj) 

i=o 


dlj{x) 

dx 


(2.2) 


(2.3) 


In  this  study,  we  will  employ  the  Chebyshev-Gauss-Lobatto  quadrature  points,  namely, 


Xi 


(2.4) 


which  are  the  roots  of  (1  —  x‘^)T^{x)  and  Tj^ix)  is  the  A^_th  degree  Chebyshev  polynomial  of  the 
first  kind.  The  Chebyshev  interpolation  polynomial  is  given  by 


_  (-iy+i(i-x^)r^(x) 

CjN‘^{x  —  Xj) 


(2.5) 


where  Cj  =  1,  j  =  1, . . . ,  A'"  —  1  and  cq  =  Cjv  =  2. 

The  form  of  the  differentiation  matrix  D  =  in  (2.3)  can  be  found  in  [7].  The  map¬ 

ping  devised  by  Kosloff  and  Tal-Ezer  [18]  will  also  be  employed  to  enhance  the  stability  of  the 
Chebyshev  collocation  scheme  [4,  8]. 


2.2  Filters 

Spectral  methods  are  highly  efficient  and  accurate,  when  the  solution  and  its  derivatives  are 
smooth.  In  the  presence  of  discontinuities,  however,  Gibbs  Phenomenon  generates  oscillations 
that  contaminate  the  solution  and  causes  the  loss  of  the  exponential  convergence.  This  is 
explained  in  spectral  space  by  the  linear  decay  rate  of  the  coefficients  a„  of  the  global  expansion: 

OO 

u(x)  =  ^a„r„(x).  (2.6) 

n=0 

In  such  a  situation,  one  can  modify  the  global  expansion  coefficients  to  enhance  the  conver¬ 
gence  properties  of  the  approximation  via  a  filter  function  (T{ri)  [26]  with  the  following  properties 

a{ri)  =  ai-r]),  cj(±l)  =  0  , 

cj(0)  =  1,  (t('')(0)=0  q=l,...,p-l 


(2.7) 
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If  a{ri)  has  at  least  p—1  continuous  derivatives,  cj(ry)  is  termed  a  p_th  order  filter.  It  was  proved 
in  [26]  that  the  hltered  expansion  converges  faster  to  the  correct  solution  than  the  unfiltered 
original  one  in  the  case  of  a  discontinuous  function.  Moreover,  the  convergence  rate  depends 
solely  on  the  order  and  the  compactness  of  the  hlter  function  (j{rj)  and  the  distance  from  the 
discontinuities. 

The  filter  function  used  in  this  study  is  the  Exponential  hlter  given  by 

cr(r/)  =  expiy—arf)  ,  (2.8) 

where  a  =  —  ln(e)  and  e  is  the  machine  zero. 

The  spectral  hltering  can  be  expressed 

•  in  the  physical  space  as; 

N 

Fryu{xi)  =  ^ji 

i=o 


•  in  the  transformed  space  as: 

N 

Fffu{xi)  =  '^ancr  Tn{xi), 

n=0 


(2.10) 


where  Fpf  is  the  hltering  operator. 

While  hltering  techniques  can  improve  the  overall  convergence  properties  of  the  solution  away 
from  discontinuities,  the  solution  near  the  discontinuities  remains  poor.  Post-processing  of  the 
resulting  oscillatory  data  to  recover  spectrally  accurate  non-oscillatory  results  can  be  performed 
by  various  reconstruction  techniques,  such  as  the  direct  and  inverse  Gegenbauer  reconstruction 
[11,  16]  and  Fade  reconstruction  [21].  Reconstruction  techniques  are,  in  general,  computationally 
costly  and  certain  complications  might  arise.  For  instance,  as  the  degree  of  the  reconstructed 
polynomial  increases,  the  Gegenbauer  transformation  matrices  become  ill-conditioned  due  to 
the  round-off  error  [16].  Moreover,  the  extension  of  these  reconstruction  techniques  to  higher 
dimensions  is  not  trivial. 

2.3  The  WENO  Method 

WENO  schemes  were  designed  for  the  numerical  solution  of  Hyperbolic  Conservation  Laws  in 
the  form  (for  ease  of  presentation  the  discussion  is  based  on  the  one-dimensional  formulation) 

ui  +  f(u),  =  0  .  (2.11) 

The  Jacobian  for  (2.11)  is  needed  to  project  the  system  into  its  characteristic  form  and, 
in  general,  ^  is  not  constant.  To  overcome  the  difficulty  when  computing  the  flux  at  a  cell 
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boundary  ’’local  freezing”  of  the  matrix  components  is  employed  using  the  Roe  average 

Uj+i  [23]  defined  as 

/(ui+i)  -  /(uj)  =  /'(u._^i)(uj+i  -  Ui)  .  (2.12) 

The  numerical  scheme  for  (2.11)  can  be  written  in  the  conservative  form  as 


dui{t) 

dt 


Ax  ^'^*+2 


(2.13) 


where  Ui{t)  is  the  numerical  approximation  to  the  point  value  u{xi,t)  and  Ax  is  the  uniform 
grid  spacing.  The  numerical  flux  based  on  the  {uj,j  =  i  —  r, . . .  +  s} 


f  {Ui—r:  •  •  •  )  lli+s) 


(2.14) 


satisfies  the  following  conditions: 

•  /  is  Lipschitz  continuous  in  all  arguments; 

•  /  is  consistent  with  the  physical  flux  /,  i.e.  f{u, . . .  ,u)  =  f{u). 

The  solution  to  the  conservative  scheme  (2.13),  if  converges,  will  converge  to  a  weak  solution. 
This  is  known  as  the  Lax-Wendroff  theorem. 

The  earlier  works  of  van  Leer  [25],  Boris  and  Book  [1]  and  Harten’s  TVD  schemes  have  led 
to  the  introduction  of  Essentially  Non-Oscillatory  (ENO)  schemes  [14],  where  its  basic  premise 
is  that  of  adaptivity  of  the  stencil,  based  on  the  local  smoothness  of  the  solution,  using  only 
a  stencil  free  of  discontinuities  in  the  computation  of  flux  gradients.  ENO  schemes  have  been 
shown  to  generate  sharp  resolutions  of  shocks  as  well  as  to  maintain  high  order  of  accuracy  in 
smooth  regions. 

An  improvement  of  finite  volume  ENO  -  Weighted  Essentially  Non-Oscillatory  (WENO) 
schemes  were  proposed  by  Liu,  Osher,  and  Chan  [19].  WENO  schemes  use  the  same  idea  as 
ENO,  except  that  WENO  uses  a  convex  combination  of  all  available  smooth  stencils  to  obtain 
higher  order  of  accuracy  than  the  original  ENO  at  smooth  parts  of  the  solution.  Finite  Difference 
WENO  schemes  were  later  proposed  by  Jiang  and  Shu  [15]. 

In  all  numerical  examples  that  follows,  we  use  the  fifth  order  characteristic-wise  WENO  finite 
difference  formulation,  making  use  of  the  Roe  average  for  the  eigensystem  defined  above  and 
the  global  Lax-Fredrichs  flux  splitting 

±au)  ,  (2.15) 

where  a  =  max^  maxi<j<Ar^  l^j('*^)l  and  Xi{u),i  =  l,...,Nx  are  the  local  eigenvalues  of  the 
Jacobian  The  numerical  flux  /j^i  is  computed  by  the  WENO  reconstruction  procedure. 
The  interested  reader  is  referred  to  [3]  for  details. 
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3  Multi-Resolution  Analysis 

The  successful  implementation  of  the  Hybrid  method  depends  on  the  ability  to  obtain  accurate 
information  on  the  smoothness  of  a  function.  In  this  work,  we  employ  the  Multi-Resolution 
(MR)  algorithms  by  Harten  [12,  13]  to  detect  the  smooth  and  rough  parts  of  the  numerical 
solution.  The  general  idea  is  to  generate  a  coarser  grid  of  averages  of  the  point  values  of  a 
function  and  measure  the  differences  (MR  coefficients)  di  between  the  interpolated  values  from 
this  sub-grid  and  the  point  values  themselves.  A  tolerance  parameter  eMR  is  chosen  in  order 
to  classify  as  smooth  those  parts  of  the  function  that  can  be  well  interpolated  by  the  averaged 
function  and  as  rough  those  where  the  differences  di  are  larger  than  the  parameter  Cma-  We 
shall  see  that  the  order  of  interpolation  is  relevant  and  the  ratio  between  di  of  distinct  orders 
may  also  be  taken  as  an  indication  of  smoothness. 

Let  us  start  by  showing  two  examples  where  one  can  notice  the  detection  capabilities  of  the 
Multi-Resolution  analysis  that  will  be  presented  below.  The  left  and  right  figures  of  Figure  1 
show  the  piecewise  analytic  function 

(  10 -|- — 1  <  X  <  — 0.5 

/(x)  =  I  x^  -0.5  <  X  <  0  ,  (3.1) 

sin(27rx)  0  <  x  <  1 

and  the  density  {p)  of  the  Mach  3  Shock-Entropy  wave  interaction  problem  [15]  as  computed  by 
the  classical  fifth  order  WENO  finite  difference  scheme,  respectively. 


X 


Figure  1:  (Left)  The  third  and  eighth  order  MR  coefficients  di  of  the  piecewise  analytic  function.  (Right) 
The  third,  fifth  and  seventh  order  MR  coefficients  di  of  the  density  /(x)  =  p  of  the  Mach  3  Shock-Entropy 
wave  interaction  problem. 


The  test  function  (3.1)  has  a  jump  discontinuity  at  x  =  —0.5  and  a  discontinuity  at  its  first 
derivative  at  x  =  0.  One  can  see  that  at  each  grid  point  the  differences  di  decay  exponentially  to 


7 


Costa  et.  al  /  Commun.  Comput.  Phys.,  1  (2006),  pp.  1-31 


zero  inside  the  analytical  pieces  of  the  function  when  the  order  of  interpolation  increases  from 
n^MR  =  3  to  Umr  =  8.  At  the  discontinuity  x  =  0.5,  the  measured  differences  di  are  0(1)  and 
remain  unchanged  despite  the  increase  of  the  interpolation  order.  Similar  behavior  is  exhibited 
at  the  derivative  discontinuity  at  x  =  0  with  a  smaller  amplitude. 

Also,  in  the  right  hgure  of  Figure  1,  the  density  of  the  Mach  3  Shock-Entropy  wave  interaction 
problem  and  the  corresponding  MR  coefficients  di  are  shown  for  the  third,  fifth  and  seventh 
order  Multi-Resolution  analysis.  The  location  of  the  main  shock  is  at  x  ~  2.73  and  the  shocklets 
behind  the  main  shock  are  well  captured.  The  high  frequencies  behind  the  main  shock  are  much 
better  distinguished  with  the  higher  orders. 

Averaging  a  function  corresponds  to  filter  the  upper  half  of  the  spectrum.  The  main  idea 
of  Hartens  smoothness  classihcation  is  to  measure  how  distant  the  actual  values  of  the  function 
are  from  being  predicted  through  interpolation  of  the  lower  half  of  the  frequencies  contained  in 
the  sub-grid  of  averages.  We  now  describe  a  detailed  construction  of  the  sub-grid  of  averages 
and  its  corresponding  interpolating  polynomial,  finishing  with  a  worked  example. 

Given  an  initial  number  of  grid  points  Nq  and  grid  spacing  Axq,  consider  the  set  of  nested 
dyadic  grids  {G^,0  <  k  <  L},  defined  as: 


=  {xti  =  0,...,Nk},  (3.2) 

where  xf  =  iAxfc,Axfc  =  2^Axo,Nk  =  2~^Nq.  For  each  level  k  >  0  we  define  the  set  of  cell 
averages  {f^,  i  =  1, . . . ,  N^}  at  x^  of  a  function  /(x): 

fi  =  [  fix)dx,  (3.3) 

Axfc  J^k 

t  —  1 


and  =  f^.  Let  f^i-i  be  the  approximation  to  by  the  unique  polynomial  of  degree  2s 

that  interpolates  |(|  <  s  at  where  r  =  2s  -|-  1  is  the  order  of  approximation. 

The  approximation  differences,  also  called  multiresolution  coefficients,  d\  =  /^jTi  — 
at  the  A:_th  grid  level  and  grid  point  Xj,  have  the  property  that  if  /(x)  has  p  —  1  continuous 
derivatives  and  a  jump  discontinuity  at  its  p_th  derivative,  then 


Ax^[/}^^]  for  p  <r 
for  p  >  r 


(3.4) 


where  [•]  denotes  the  magnitude  of  the  jump  of  the  function  inside. 

From  formula  (3.4)  it  follows  that 

~  I,  where  p  =  min{p,r}.  (3.5) 

Equation  (3.5)  shows  that  away  from  discontinuities,  the  MR  coefficients  d^  diminish  in  size 
with  the  refinement  of  the  grid;  close  to  discontinuities,  they  remain  the  same  size,  independent 
of  k.  The  MR  coefficients  d^  were  used  in  [13]  in  two  ways.  First,  finer  grid  data  were  mapped 
to  its  M  level  multiresolution  representation  ff  =  {d},  -  ■  ■  ,df ,  fl^)  to  form  a  multiscale  version 
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of  a  particular  scheme,  where  truncation  of  small  quantities  with  respect  to  a  tolerance  parameter 
decreased  the  number  of  flux  computations.  Secondly,  the  MR  coefficients  also  acted  as  a 
shock  detection  mechanism  and  an  adaptive  method  was  designed  where  a  second-order  Lax- 
Wendroff  scheme  was  locally  switched  to  a  first-order  accurate  TVD  Roe  scheme,  whenever  dj 
was  bigger  than  eMR- 

Equation  (3.4)  also  indicates  that  the  variation  of  the  MR  order,  Umr,  can  give  additional 
information  on  the  type  of  the  discontinuity.  Nevertheless,  in  this  work,  we  will  be  limited 
at  using  only  the  first  level  A:  =  1  of  the  multiresolution  coefficients  and  we  shall  drop  the 
superscript  1  from  the  dj  from  here  on  unless  noted  otherwise. 

Hence,  to  find  di,  the  idea  is  to  construct  a  piecewise  polynomial  Pk{x)  of  degree  k  =  Umr 
using  k  +  1  computed  average  values  of  fi,  fi,  at  the  equi-spaced  grid  Xi  such  that 

Pkixi)  =  f{xi)  +  0(Ax^+^),  (3.6) 

and 

di  =  fi- Pk{xi).  (3.7) 

Given  a  tolerance  level  Cmh,  the  smoothness  of  the  function  f{x)  at  Xi  would  then  be  checked 
against  the  magnitude  of  the  di,  namely: 

/  \di\  <  eMR  solution  is  smooth.  ,  . 

\  \di\  >  Cm R  solution  is  non-smooth.  ^  ' 

The  algorithm  for  computing  the  MR  coefficients  di  is  given  next. 

3.1  Computing  the  MR  Coefficients 

Consider  an  equi-spaced  grid  {xi  =  zAx,  i  =  — m, . . .  ,0, . . . ,  N, . . . ,  N  +  M}  where  Ax  is  the 
constant  grid  spacing.  N  can  be  an  odd  or  an  even  number.  Depending  on  N  and  the  even  or 
odd  order  of  the  MR  Analysis  Umr,  the  number  of  ghost  points  m  and  M  required  are  given  in 
table  I. 


N 

klMR 

m 

M 

odd 

odd 

XI-MR  +  1 

TXmR  +  1 

even 

odd 

^mr  +  1 

TXmR 

odd 

even 

klMR 

eiMR  +  2 

even 

even 

klMR 

eiMR  +  1 

Table  I:  The  number  of  ghost  points  m  and  M  required  for  the  MR  Analysis. 


Given  the  grid  point  values  of  the  function  /(x),  the  average  values  are  computed  as 
J  ^  if  ,  f  ^  N  +  M  -1 

P  —  2  ^  “  — 2'  "  '  ’ - ' 


2 


(3.9) 
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We  construct  a  piecewise  k  =  Umr  degree  polynomial  Pk{x)  using  the  k  +  1  computed  average 
values  of  the  given  function,  /*  such  that 

Pk{xi)  =  f{xi)  +  0(Ax^+^).  (3.10) 


The  polynomial  Pk{xi),  I 
be  written  as 


and  L  =  (  —  1  or  L  =  (  if  /c  is  odd  or  even,  respectively,  can 
i-\-L 

Pk{Xi)  =  ^  arjr-  (3.11) 

r=i—l 


However,  since  the  coefficients  a  depend  only  on  Xi  and  do  not  depend  on  the  function  /(x), 
the  Pk{xi)  can  be  written  as: 


Pdx.)  =  l 

[  2^r=-l 

mod(i,  2)  =  0 
mod(i,  2)  =  1 

(3.12) 

In  the  case  of  mod(i,  2)  =  1, 

(3-r  =  ar,  r  =  —1, 

...,L. 

(3.13) 

Furthermore,  if  Umr  is  even,  the  coefficients  a  are  symmetric  about  r  =  0,  namely,  a-r 

II 

s- 

e 

II 

1, . . . ,  L. 

The  desired  coefficients  a  are  computed  by  requiring  Pk{x)  to  be  equal  to  each  of  the  first 
k  +  1  monomials  /(x)  =  1,  x,  x^, . . . ,  x^  and  evaluated  at  any  grid  point  x  =  x* .  For  simplicity, 
we  take  x*  =  0.  The  /*  are  evaluated  for  i  =  —I, . . . ,  L.  This  procedure  results  in  a  system  of 
linear  equations,  Acf  =  b,  where 


A  = 

/  1 

-21  +  {-21  -hi) 

1  \ 
2L  -h  (2L  -h  1) 

,  a  = 

/  a_i  \ 

a-z+i 

,~b  = 

(  1  \ 
0 

(-20^  +  (-2(  +  1)^  ., 

. .  (2L)^  +  (2L  +  1)^  yi 

V  ttL  / 

V  0  / 

(3.14) 


and  A  is  a  matrix  of  size  (L  1)  x  (L  1). 

Using  (3.12),  the  /c-th  order  Multi-Resolution  coefficients  d*  at  x*  can  be  computed  as 


di  =  fi-  Pk{xi)  i  =  0, . . . ,  A. 


(3.15) 


One  can  also  evaluate  the  a  by  matching  the  terms  in  the  Taylor  series  expansion  using 
(3.10)  and  (3.11)  to  any  desired  order,  however  this  procedure  may  become  cumbersome  for 
high  order  k. 


Example 


To  illustrate  the  procedure  above,  we  will  construct  two  unique  local  polynomials  with  k  = 
nuR  =  3,  such  that  Pkixo)  =  /(xq)  -|-  0{Ax^~^^)  and  Pk{xi)  =  f{xi)  -|-  0{Ax^~^^). 
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To  construct  the  desired  polynomials  one  needs  to  find  the  unique  coefficients  {a-2,  ct-i,  cto,  cti 
and  {/3_i, /3o, /?!, /32}  such  that 


a_2/-2  +  a-i/-i  +  ao/o  +  ai/i  =  f{xo)  +  0(Ax'‘) 


and 


P-if-i  +  Pofo  +  Pifi  +  P2f2  =  /(ail)  +  O(Ax^). 

The  system  of  equations,  (3.14),  becomes 


A  = 


/I  1  1  1  \ 

/  a-2  \ 

-7-315 

a_i 

~h  — 

0 

25  5  1  13 

5  Of  - 

ao 

?  0  — 

0 

\  -91  -9  1  35  ) 

\  «i  / 

VO/ 

Solving  this  system  yields 


17 


55 


(3.16) 

(3.17) 


(3.18) 


(3.19) 


and  {P-i  =  ai,Po  =  ao,Pi  =  a-i,P2  =  a-2}- 


Remark  3.1.  The  tolerance  parameter  Emr  determines  the  dynamic  activation  of  the  spectral 
and  WENO  spatial  discretizations  along  the  various  subdomains  of  the  hybrid  method.  While 
a  too  small  value  of  Cmr  activates  the  more  expensive  WENO  method  at  subdomains  where  the 
solution  is  smooth,  a  larger  value  activates  the  spectral  method  at  a  subdomain  with  low  spatial 
resolution,  generating  oscillations.  €mr  also  bears  a  straight  relation  with  the  interpolation  order 
Umr-  High  Umr  values  decrease  the  size  of  Emr  one  needs  to  chose,  since  high  frequencies  are 
less  mistaken  by  gradient  jumps.  The  general  guideline  is  to  start  with  a  value  for  Umr  at  least 
equal  to  the  order  of  the  WENO  method  and  increase  it  according  to  the  complexity  of  the 
solution.  Eor  instance,  Umr  =  5  is  a  good  choice  for  the  piecewise  smooth  solution  of  the  SOD 
problem,  the  Entropy  problem  would  work  better  with  Umr  =  7.  For  most  of  the  flows  with 
shock  that  were  tested,  the  value  of  Cma  =  10“^  yielded  a  good  balance  between  computational 
speed  and  accuracy  of  the  numerical  solution. 


4  The  Multi-Domain  Hybrid  Spectral- WENO  method 

We  now  describe  the  implementation  of  the  Multi-Domain  Hybrid  Spectral- WENO  method  (Hy¬ 
brid),  detailing  the  structure  of  the  two-dimensional  grid  of  subdomains,  the  interfaces  treatment 
and  the  algorithm  that  switches  the  spatial  discretization  of  the  subdomains. 

The  main  idea  of  the  Hybrid  method  can  be  formulated  as: 

Avoid  Gibbs  phenomenon  by  keeping  diseontinuities  at  WENO  subdomains  and  inerease  the 
numerical  efficiency  by  treating  the  smooth  parts  of  the  solution  using  a  spectral  spatial  dis¬ 
cretization. 
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Figure  2:  (Left)  Partition  of  the  physical  domain  into  spectral  and  WENO  subdomains.  (Middle)  Typical 
spectral  subdomain.  (Right)  Typical  WENO  subdomain. 

In  this  study,  the  physical  domain  is  restricted  to  rectangular  shapes  and  will  be  partitioned 
into  a  {N^  x  grid  of  subdomains.  Figure  2  shows  an  example  of  such  a  domain  partition 
along  with  typical  spectral  and  WENO  subdomains.  Note  that  patching  of  subdomains  occurs 
due  to  the  ghost  points  of  the  WENO  discretization. 

We  shall  use  a  vector  k  =  {kx,  ky),kx  =  1,  •  •  • ,  ky  =  1, . . . ,  Wj  to  denote  the  coordinates 
of  the  two  dimensional  subdomain  grid.  For  example,  k  =  (2, 3)  means  the  subdomain  number 
kx  =  2  in  the  x  direction  and  =  3  in  the  y  direction.  Each  subdomain  is  initialized  either  as  a 
spectral  or  WENO  subdomain.  The  Chebyshev-Gauss-Lobatto  points  are  used  for  the  spectral 
discretization  and  an  uniformly  spaced  grid  with  ghost  points  is  used  as  the  WENO  grid. 

For  the  sake  of  simplicity,  we  consider  only  square  subdomains,  with  Ns  x  Ns  grid  points  at 
a  pectral  subdomain  and  N^r  x  N^r  equi-spaced  points  at  a  WENO  subdomain.  The  number 
of  ghost  points  is  denoted  by  r  and  Ng  is  the  number  of  points  used  for  the  buffer  zone  (see 
below) . 

4.1  Description  of  WENO  subdomains 

Each  WENO  subdomain  k  is  composed  of  three  parts;  The  ’’Ghost  Area”,  the  ’’Buffer  Area”, 
and  the  ’’Interior  Area”  as  shown  in  Figure  3. 

•  Ghost  Area: 

The  ’’Ghost  Area”  is  used  for  the  WENO  Reconstruction  and  for  communication  with  its 
neighboring  subdomains  and  is  subdivided  into  eight  ’’Ghost  Zones” 

{Gf,  i  =  1, . . . ,  8}  (see  the  middle  figure  of  Eigure  3), 

r  =  4"  X  J|  ,  =  If  X  J|  ,  Gl  =  If  X  J|  ] 

{  Gf  =  4xJf  ,  ,  G|  =  /fxJf  \ 

[g^  =  4xJ^  ,  G^  =  lfxJ^  ,  G^,=llxJl} 


(4.1) 
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Figure  3:  (Left)  Typical  WENO  subdomain.  (Middle)  ’’Ghost  Area”.  (Right)  ’’Buffer  Area”. 


where  the  stencil  ranges  are 

1q  {^0  ?  •  •  •  )  ?  ^2  {^W-t-1  ?  *  *  *  5  ^  N+r  }  : 

Jq  =  {y-r,---,y-i},  Jf  =  {yo,---,yN},  Jl  =  {yN+i,- ■  ■  ,yN+r}, 
with  N  =  Nw 


•  Buffer  Area: 

The  Buffer  Area  is  used  for  the  treatment  of  moving  discontinuities.  If  a  high  gradient 
or  shock  is  detected  inside  a  buffer  zone,  the  closest  neighboring  spectral  subdomain(s)  is 
switched  to  WENO  subdomain(s). 

For  two-dimensional  problems,  each  WENO  subdomain  has  eight  ’’Buffer  Zones” 

{Bf^,  i  =  A, ,  H}  (see  the  right  hgure  of  Figure  3), 


=  =  ,  B^  =  IIx4  ] 

Bl  =  I^oXJi  ,  ,  Bl  =  llxj\  ^ 

Bl  =  llxJl  ,  B^  =  I\x4  ,  Bl=llxJl  ] 


where  the  stencil  ranges  are  given  by 


{Xq,  .  .  .  ,  Xm};  I\  —  {Xm+1^  ■  ■  ■  )  S^AT-m})  I2  ~  {.Xn-M+Ii 
{2/0)  ■  •  ■  )  I/m})  J\  =  {Vm+Ii  ■  ■  ■  1  I/jV-m})  J2  ~  {Vn-M+Ii 


■  ■  )  Xn}, 

■  ■  ,2/iv}, 


(4.2) 


with  N  =  Niy  and  M  =  Ng.  As  greater  is  the  value  of  Ng,  earlier  is  the  detection  of  shocks 
and  gradients,  however,  at  the  cost  of  early  switching  of  the  neighboring  subdomains  to 
WENO  and  greater  chance  of  unnecessary  costly  computations.  Satisfactory  results  have 
been  obtained  with  the  default  M  =  r  in  this  study.  It  should  be  noted  that  these  buffer 
zone  grid  points  are  part  of  the  interior  grid  points  and  should  not  be  confused  with  WENO 
ghost  points. 
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•  Interior  Area: 

Finally,  the  ’’Interior  Area”  is  all  the  WENO  grid  points  excluding  the  ’’Ghost  Area” 
and  ’’Buffer  Area”,  =  W^/ (IJj  G\  0  Uj  gradient  stays  inside  this  area, 

no  action  is  required  by  the  neighbors. 

4.2  Interface  Conditions 

The  following  conhgurations  are  representative  of  any  two-dimensional  domain  partitions: 

•  Spectral-Spectral-Spectral-Spectral  (Figure  4) 

•  Spectral-Spectral-Spectral- WENO  (Figure  6) 

•  Spectral-Spectral- WENO- WENO  (Figure  6) 

•  Spectral- WENO- WENO- WENO  (Eigure  6) 

•  WENO- WENO- WENO- WENO  (Eigure  5) 

The  functional  values  in  the  spectral  and  WENO  subdomains  are  denoted  by  sfj  and  wfj 
respectively  at  {xi,yj)  in  the  respective  subdomain  k.  Centered  at  the  subdomain  k,  we  dehne 
the  superscripts  ko ,  ki ,  k2 ,  ks ,  k4 ,  ks ,  kg ,  ky ,  kg  as 


k  -|-  (—1,  4-1) 

,  ky  —  k  4-  (4-0,  4-1) 

,  kg  =  k  4- (4-1,  4-1)  \ 

k  -|-  (—1,  +0) 

,  kg  =  k 

,  kg  =  k  4-  (4-1,  4-0)  1  , 

(4.3) 

k  -|-  (—1,  —1) 

,  kg  =  k  4-  (4-0,  —1) 

,  k4  =  k  4- (4-1,  —  1)  / 

and  they  denote  the  Center  (kg).  Left  (ki),  Bottom-Left  (k2).  Bottom  (kg),  Bottom-Right  (k4). 
Right  (kg),  Top-Right  (kg).  Top  (ky),  and  Top-Left (  kg)  subdomains  of  subdomain  k  =  kg 
respectively.  Eor  example  is  to  be  interpreted  as  a  value  at  the  point  {xi,yj)  in  the  spectral 
Left  subdomain  ki  =  k  -|-  (— 1,-|-0).  The  subdomain  interfaces  consist  of  corners  and  shared 
sides  among  the  spectral  and  WENO  subdomains. 

Spectral-Spectral  Interface 

Corners  only  appear  at  the  connection  of  two  or  more  spectral  subdomains,  since  the  grid  points 
of  the  WENO  subdomains  never  coincide  with  the  spectral  collocation  points  at  the  interface. 
Corner  values  are  assigned  with  the  average  of  the  values  of  all  connecting  subdomains.  Eor 
example,  in  figure  4,  assuming  that  is  the  reference  subdomain  kg  and  N  =  Ng, 

-1-  -l-  -l-  ^  (4  4) 

•’OJV  •’OO  ^NO  ^NN  ^  y^ON  W  Oqq  “T  Ojyg  “T  Ojyjy  ^ 

Along  the  interface  between  two  spectral  subdomains,  the  values  at  the  shared  side  are 
computed  similarly  by  assigning  the  average  of  the  two  values  of  the  connected  subdomains  at 


Costa  et.  al  /  Commun.  Comput.  Phys.,  1  (2006),  pp.  1-31 


14 


1 

s 

s’ 

Figure  4:  (Left)  Spectral-Spectral-Spectral-Spectral  subdomains  configuration.  (Middle)  Corner  point  of 
all  four  spectral  subdomains.  (Right)  Shared  side  of  two  spectral  subdomains. 

all  collocation  points.  Consider  the  interface  between  the  two  spectral  subdomains  —  S'^  (the 
rightmost  figure  of  Figure  4), 

So°  =  si]  =  ^  ’  J  =  0, . . . ,  iV  .  (4.5) 

The  other  spectral  subdomains’  corners  and  sides  are  treated  similarly. 

Note  4.1.  In  case  of  adjacent  subdomains  with  different  numbers  of  collocation  points,  global 
interpolation  can  be  used  before  averaging. 

Note  4.2.  An  exact  Riemannn  solver  and  penalty  interface  conditions  had  been  tried  but  no 
discernible  differences  were  observed  from  the  simple  average. 


WENO-WENO  Interface 


To  maintain  conservation,  adjacent  WENO  subdomains  are  required  to  have  the  same  spacing 
Ax  (see  [24]).  This  implies  that  ghost  points  of  a  WENO  subdomain  match  the  interior  points  of 
the  neighboring  WENO  subdomain  (see  Figure  5).  Thus,  simple  copying  of  the  solution  values 
of  the  neighboring  interior  points  to  the  ’’Ghost  Area”  is  sufficient  for  the  communication  of 
WENO  subdomains. 

For  instance  the  middle  figure  of  Figure  5  shows  two  ’’Ghost  Zones”:  The  zone  in  the  solid 
rectangle  corresponds  to  the  ’’Ghost  Zone”  G\  of  and  the  zone  in  the  dash-dotted  rectangle 
corresponds  to  the  ’’Ghost  Zone”  of  (see  4.1).  In  this  case,  by  denoting  N  =  and 
using  as  a  reference  subdomain  (subdomain  k  =  kp) 


w 


ko 

(N+l)j 


= 


w 


ks 

-ij 


=  w 


ko 

(jV+l-i)j 


j  =  -r,. . .  ,N  +  r, 


1  =  1,. ..,r 


(4.6) 


where  r  is  the  number  of  ghost  points. 


Note  4.3.  In  general  the  number  of  ’’ghost”  points  in  each  subdomain  does  not  have  to  be  the 
same  as  long  as  the  same  grid  spacing  Ax  is  used  at  all  the  neighboring  subdomains. 
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Figure  5:  (Left)  WENO-WENO-WENO-WENO  subdomain  configuration.  (Middle)  Interface. 

(Right)  Interface. 


Spectral- WENO  Interface 


Spectral-WENO  interface  can  be  found  in  the  configurations  shown  in  Figure  6. 


Figure  6:  Subdomain  Configurations:  (Left)  Spectral-Spectral-Spectral- WENO.  (Middle)  Spectral- 
Spectral- WENO- WENO.  (Right)  Spectral- WENO- WENO- WENO. 


The  ’’Ghost  Zones”  and  of  the  WENO  subdomain  lie  in  the  three  neighboring 

spectral  subdomains  as  shown  in  Figure  6.  The  ’’Ghost  Area”  points  of  WENO  subdomains 
are  computed  via  spectral  interpolations.  To  be  more  specific,  referring  to  the  middle  figure  of 
Figure  7,  the  solution  in  the  area  of  the  intersection  of  (subdomain  k)  and  (subdomain  ks) 
corresponding  to  the  ’’Ghost  Zone”  G^  in  (4.1)  (see  also  Figure  3)  is  read  using  the  interpolating 
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Figure?:  (Left)  Spectral-Spectral-Spectral- WENO  subdomain  configuration.  (Middle)  5^  — IE'*  Interface. 
(Right)  Interface. 

polynomial  representing  the  solution  in  S^,  that  is,  at  the  subdomain  k, 

Ns  Ns 

'wiXn,ym)  =  '^'^Sijli{Xn)lj{ym),  m  =  0,...,Nw,  n  =  -r,  .  .  .  , -1,  (Xn,t/m)Gk5  , 
j=0  i=0 

(4.7) 

where  Sij  are  the  functional  values  in  (subdomain  k),  li(x)  and  lj{y)  are  Lagrangian  inter¬ 
polation  polynomials  of  degree  Ns-  The  ghost  points  for  the  other  two  configurations  can  be 
obtained  in  a  similar  way. 

The  interface  points  of  the  spectral  subdomain  are  computed  via  a  two  dimensional  interpola¬ 
tion  polynomial  of  degree  r  =  n^.  The  choice  of  stencils  {{xi,yj),i  =  —iQ  —  rf2, . . .  ,iQ  +  r/2,j  = 
jo  —  r/2, . . .  ,jo  +  r/2}  for  the  interpolation  polynomial  should  be  as  symmetric  about  a  given 
point  ixiQ,yj^)  as  possible. 

4.3  The  Switching  Algorithm 

In  this  section  we  describe  the  algorithm  used  to  switch  the  subdomains  spatial  discretizations 
as  indicated  by  the  Multi-Resolution  Analysis  of  Section  3.  The  following  three  conditions  are 
the  main  rules  to  be  followed: 

1.  If  a  subdomain  contains  high  gradients,  then  switch  its  spatial  discretization  to  (or  keep 
it  with)  WENO; 

2.  If  high  gradients  are  present  in  the  ’’Buffer  Areas”  of  connected  neighboring  subdomains, 
then  switch  the  current  subdomain  to  (or  keep  it  as)  a  WENO  subdomain; 

3.  In  any  other  case,  switch  the  subdomain  to  (or  keep  it  as)  a  spectral  subdomain; 

The  first  condition  above  avoids  the  Gibbs  phenomenon,  keeping  the  discontinuities  inside 
WENO  subdomains.  The  second  condition  ensures  the  switch  to  WENO  subdomain  in  order  to 
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allow  only  WENO-to-WENO  transmission  of  discontinuities.  The  third  condition  improves  the 
numerical  efficiency,  since  it  ensures  that  smooth  parts  of  the  solution  will  always  be  contained 
in  spectral  subdomains. 

Multi-Resolution  analysis  of  the  solution  is  performed  at  every  sub-stage  of  the  third  order 
TVD  Runge-Kutta  scheme  used  for  the  temporal  evolution.  At  each  subdomain  k,  we  dehne 
the  smoothness  flag  variable,  FlagJ^-,  at  each  grid  point  {xi,yj)  including  the  ghost  points,  as  , 


1 1.  m  >  ^MR  for  {xi,yj)  C  (  =  A,...,H 

\  0,  otherwise 


(4.8) 


where  dfj  are  the  MR  coefficients.  Since  the  Multi- Resolution  Analysis  requires  uniformly  spaced 
grids,  the  spectral  grids  are  hrst  interpolated  to  uniformly  spaced  grids  before  obtaining  the  MR 
coefficients  dfj.  The  necessary  ghost  points  are  acquired  from  the  neighboring  subdomains.  At 
the  boundary  subdomains,  the  value  of  the  boundary  ghost  points  are  extrapolated  linearly  from 
the  interior  data. 

The  algorithm  proceeds  by  checking  for  each  subdomain  k  and  at  the  buffer  zones  of  the 
neighboring  subdomains  if  Flag^-  is  equal  to  one.  If  so,  it  switches  subdomain  k  to  a  WENO 
discretization.  On  the  other  hand,  if  FlagJ^-  is  identically  zero,  then  a  spectral  discretization  is 
implemented,  or  kept.  These  switches  require  the  use  of  interpolation  from  a  Chebyshev  grid  of 
points  to  a  uniformly  spaced  one  and  vice-versa: 


•  To  switch  from  the  spectral  subdomain  to  the  WENO  subdomain,  the  data  are  interpolated 
onto  the  uniformly  spaced  grid  via  the  spectral  interpolation  formula. 


•  To  switch  from  the  WENO  subdomain  to  the  spectral  subdomain,  the  data  are  interpolated 
onto  the  Chebyshev  Gauss-Lobatto  points  via  the  Lagrangian  interpolation  polynomial  of 
the  same  order  as  the  WENO  method. 


Remark  4.1.  Back  and  forth  switching  between  WENO  and  spectral  discretizations  may  occur 
too  frequently  for  the  same  domain  when  the  Cma  is  marginally  set.  The  dfj  coefficients  might 
oscillate  around  the  parameter  €mr  in  time  due  to  some  numerical  factors  such  as  dissipation, 
dispersion  and  nonlinear  effects,  or  any  combination  of  such.  This  pattern  of  switching  can 
repeat  itself  for  a  while  until  the  solution  settles  down  with  a  clear  dehnition  of  the  dfj,  which  is 
either  greater  than  or  smaller  than  the  MR  tolerance  Cmr-  In  order  to  alleviate  such  occurrences, 
one  must  devise  a  procedure  preventing  the  switch  from  WENO  to  spectral  if  it  had  already 
occurred  recently.  However,  the  procedure  must  never  prevent  a  spectral  to  WENO  switch,  for 
oscillations  and  instability  might  occur. 


5  Numerical  Results 

In  this  section,  we  apply  the  hybrid  method  to  two  well-known  problems  in  Conservation  Laws; 
The  Shock- Vortex  Interaction  and  the  Richtmyer- Meshkov  Instability.  The  governing  equations 
are  the  two-dimensional  Euler  equations  in  Cartesian  coordinates  given  by: 

Qi  +  Fa,  -|-  Gy  =  0 


(5.1) 
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where 


Q  =  {p,pu,pv,Ef  , 

F  =  {pu,pu^  +  P,puv,{E  +  P)u)'^  ,  (5.2) 

G  =  {pv,  puv,  pv^  +P,{E  +  P)vY  , 

and  the  Equation  of  state 

P  =  il +  ]^p{Y +  vY),  7  =  1-4.  (5.3) 

We  will  restrict  the  following  study  to  rectangular  domains,  which  will  be  partitioned  into  an 
equal  number  of  subdomains  in  both  x  and  y  directions.  The  number  of  Chebyshev  collocation 
points  for  all  spectral  subdomains  will  be  the  same  in  both  x  and  y  directions,  as  well  as  the 
number  of  uniformly  spaced  grid  points  for  the  WENO  subdomains.  The  order  of  the  Multi- 
Resolution  Analysis  is  the  same  as  the  WENO  scheme,  Umr  =  Hw  A  14_th  order  Exponential 
filter  and  Kosloff-Tal-Ezer  mapping  are  employed  in  all  spectral  subdomains.  Eree  stream  bound¬ 
ary  conditions  are  imposed  in  the  inflow  and  outflow  in  the  x  direction  and  periodical  boundary 
condition  is  imposed  in  the  y  direction.  To  evolve  the  ODE  from  the  semi-discretized  PDE  in 
time,  the  third  order  Total  Variation  Diminishing  Runge-Kutta  scheme  (RK-TVD)  will  be  used 
[3]: 

=  U^  +  AtL{UY 

t72  =  ^  (su^  +  +  AtL{uY)  ,  (5.4) 

C7”+^  =  ^  +  2AtL{UY) 

where  L  is  the  spatial  operator.  The  CEL  numbers  for  the  spectral  and  WENO  subdomains  are 

set  to  be  3  and  0.4,  respectively.  All  numerical  experiments  were  run  on  a  667  MHz  Compaq 

Alpha  machine  with  1GB  memory  and  with  an  Alpha  internal  floating  point  processor. 

5.1  Shock  —  Vortex  Interaction 

The  tangential  velocity  profile  of  the  counter-clockwise  rotating  vortex  [17]  centered  at  {xc,yc) 
and  strength  T  is  given  in  polar  coordinates  by: 

f  Tr{rY-rY) 

U{r)  =  <  rr(r“2  —  r^) 

I  0 

where  rg  =  0.2  and  ri  =  1.0. 

Due  to  the  strong  nonlinearity  of  the  shock-vortex  interaction  its  physics  are  not  well  un¬ 
derstood.  A  better  understanding  will  have  many  potential  applications,  for  instance,  subsonic 
and  supersonic  jet  nozzle  design  and  blade  noise  generation.  The  Hybrid  method  will  be  tested 
on  this  problem  with  a  vortex  of  amplitude  T  =  0.25  and  shock  Mach  numbers  Mg  =  1.25,3,6. 


0  <  r  <  ro  <  ri 

ro  <  r  <  ri  ,  (5-5) 

r  >  ri 
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Mach  1.25 

In  this  first  example,  the  physical  domain  (0  <  x  <  3.9,  —2  <  y  <  2)  is  partitioned  into  a  13  x  10 
grid  of  subdomains.  Spectral  subdomains  use  a  32  x  32  grid  of  Chebyshev  points  and  WENO 
grids  are  50  x  50.  MR  analysis  is  performed  with  the  MR  tolerance  set  to  Emr  =  5  x  10“^. 


X 


(a)  t  =  0.73 


(b)  t  =  1.08 


(c)  t  =  1.20 


Figure  8:  Density  p  of  the  Shock- Vortex  interaction  with  Mach  number  Mg  =  1.25  and  F  =  0.25  at  (a) 
t  =  0.73,  (b)  t  =  1.08  and  (c)  t  =  1.2,  as  computed  by  the  Hybrid  method. 

When  an  initially  planar  shock  wave  hits  the  vortex,  it  deforms  as  compression  and  rarefac¬ 
tion  regions  are  created  behind  the  shock,  as  shown  in  Figure  8  for  t  =  0.732, 1.08, 1.2.  As  the 
interaction  proceeds  over  time,  strong  bifurcation  and  deformation  of  the  shock  are  observed. 
The  shock  emanating  from  the  compression  region  has  one  part  moving  upward  and  another 
moving  downward.  The  strength  of  the  lower  upward  moving  part  is  greater,  due  to  the  direc¬ 
tion  of  rotation  of  the  vortex.  MR  Analysis  performed  well  in  capturing  the  shock  and  the  high 
gradient  regions  immediately  behind  the  shock,  as  indicated  by  the  WENO  subdomains  enclosed 
with  black  bounding  boxes.  The  remaining  subdomains  are  accurately  dealt  with  by  the  spectral 
methods.  The  number  of  WENO  subdomains  is  far  fewer  than  spectral  ones  resulting  in  a  more 
efficient  algorithm  than  the  classical  WENO  scheme. 

Mach  3 

We  now  increase  the  Mach  number  of  the  shock  from  Mg  =  1.25  to  Mg  =  3  and  compare  the 
solution  of  the  Hybrid  scheme  with  a  highly  resolved  one  computed  with  the  classical  fifth  order 
WENO  scheme  using  1200  x  1200  points.  The  physical  domain  (0  <  x  <  3.0,  —2  <  y  <  2)  is 
partitioned  into  a  10  x  10  grid  of  subdomains  and  all  other  parameters  are  as  in  Example  1. 

As  in  the  previous  example,  an  acoustic  wavefront  is  generated  and  a  number  of  fine  scale 
structures  are  formed  behind  the  main  shock.  This  example  indicates  that  the  Hybrid  method 
captures  the  fine  features  of  the  solution  even  in  the  case  of  strong  shocks. 
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Mach  6 

In  this  last  example,  we  further  increase  the  shock  Mach  number  to  Mg  =  6  and  use  a  19  x  20 
grid  of  subdomains  to  partition  the  physical  domain  (0  <  x  <  4.2,  —3  <  y  <  3).  The  spectral 
grid  is  16  x  16.  This  example  shows  that  the  Hybrid  method  can  be  applied  to  higher  Mach 
number  flows  as  well,  still  resulting  on  a  reliable  and  efficient  shock-capturing  method. 


5.2  Two-dimensional  Richtmyer-Meshkov  Instability 

Richtmyer  in  1960  [22]  theoretically  predicted  the  occurrence  of  instability  on  a  perturbed  mate¬ 
rial  interface  under  the  impulsive  acceleration  of  an  incident  shock  wave.  In  1970  Meshkov  [20] 
experimentally  confirmed  these  predictions.  A  variety  of  motions  can  be  generated  following  the 
interaction  of  a  shock  wave  with  an  interface  separating  two  materials.  Any  small  perturbation 
present  on  the  interface  will  be  amplified  after  such  a  contact.  This  class  of  problems  is  referred 
in  the  literature  as  the  ’’Richtmyer-Meshkov  Instability  (RMI)”.  As  the  interface  between  two 
materials  becomes  increasingly  distorted  other  instabilities  such  as  the  Kelvin-Helmholtz  Insta¬ 
bilities  develop  and  a  region  of  turbulence  mixing  ultimately  results.  The  RMI  arises  in  many 
applications  as,  for  instance,  the  Inertial  Confinement  Fusion  (ICF)  process.  A  recent  model 
under  extensive  study  consists  of  a  set  of  laser  beams  directed  into  a  chamber  containing  a  spher¬ 
ical  fusion  fuel  target.  The  expected  result  should  be  compression,  ignition  and  a  subsequent 
energy  surplus.  However,  since  no  perfect  capsule  exists,  irregularities  on  the  surface  excite  the 
undesirable  RMI,  reducing  the  effective  uniform  compression  pressure  onto  the  capsule. 

Presently,  a  rectangular  domain  with  a  shock  Mach  number  Mg  interacting  with  a  single 
mode  sinusoidal  perturbation  along  a  Xenon  (Xe)  and  Argon  (Ar)  gases  interface  is  simulated 
using  the  Hybrid  method.  The  initial  condition  is  given  in  Figure  II  and  a  diffusive  interface  is 
modeled  with  an  exponential  function,  i.e. 


(  1  d<0 

S{x,y)  =  <  exp(— a|(i|^)  0  <  d  <  1 

1  0  d  >  1 


(5.6) 


where 


d 


{xg  +  acos(27ry/A)  -|-  d)  —  x 
2d 


(5.7) 


d  =  0.2  cm  >  0  is  the  interface  thickness,  /?  =  8  is  the  interface  order,  Xg  =  0.5  cm  is  the 
location  of  the  interface  and  a  =  —  In  e,  where  e  is  the  machine  zero.  The  conservative  or 
primitive  variables  are  scaled  according  to  the  S{x,y)  between  the  Xenon  and  Argon  gases. 


Mach  4.46 

In  this  example,  the  physical  domain  (0  <  x  <  5.1,0  <  y  <  A)  is  partitioned  into  a  17  x  12 
grid  of  subdomains.  The  spectral  grids  are  32  x  32  and  the  WENO  grids  are  50  x  50.  The  MR 
tolerance  is  now  lowered  to  €mr  =  5  x  10“^. 
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As  the  shock  wave  collides  with  the  interface  separating  the  two  gases,  the  sine  wave  pertur¬ 
bation  is  accelerated,  compressed  and  amplified  following  the  refraction  of  the  shock  (see  Figure 
13). 

The  heavier  Xenon  gas  (Xe)  will  penetrate  into  the  lighter  Argon  gas  (Ar)  forming  finger-like 
structures  -  bubbles  and  spikes.  A  bubble  (spike)  is  a  portion  of  the  light  (heavy)  gas  penetrating 
into  the  heavy  (light)  gas.  The  basic  mechanism  of  these  instabilities  is  the  baroclinic  generation 
of  vorticity  cJ  induced  from  the  misalignment  of  the  pressure  gradient  Vp  of  the  shock  and  the 
local  density  gradient  Vp  across  the  interface: 

doj 

—  ~  Vp  X  V/9,  0  =  V  X  u,  (5.8) 

where  u  is  the  velocity.  Once  again,  as  seen  in  Figure  13,  the  WENO  method  is  activated  only 
along  the  material  interface  and  where  high  gradients  appear. 

Mach  8 

The  physical  domain  (0  <  x  <  24.6, 0  <  y  <  A)  is  partitioned  into  a  82  x  12  grid  of  subdomains 
in  order  to  apply  the  Hybrid  method  for  a  longer  time  integration  up  to  t  =  lOOys.  The  shock 
Mach  number  is  increased  from  Mg  =  4.46  to  Mg  =  8  and  the  spectral  grid  is  set  to  24  x  24. 

It  can  be  observed  in  Figure  14  that  the  Hybrid  method  successfully  tracks  shocks  and 
high  gradients  with  WENO  discretization  (black  bounding  boxes)  while  the  smooth  parts  of  the 
solution  are  well  represented  using  spectral  subdomains. 

5.3  CPU  Timing 

The  Hybrid  method  has  the  potential  advantage  of  being  faster  than  the  classical  WENO  method 
due  to  the  higher  numerical  efficiency  of  spectral  methods  at  smooth  parts  of  the  solution.  More¬ 
over,  spectral  discretizations  avoid  the  expensive  characteristic  decompositions  and  projections 
of  the  WENO  method.  In  this  section,  we  provide  some  CPU  timing  results  for  the  Mach 
3  Shock- Vortex  Interaction  (see  section  5.1)  when  using  the  Hybrid  and  the  classical  WENO 
methods  with  equal  resolution. 

The  physical  domain  is  partitioned  into  a  set  of  subdomains  of  sizes  ((10  x  2-^ )  x  (10  x  2-^ )),  j  = 
0, 1,2,3.  Spectral  subdomains  use  a  12  x  12  grid  of  Chebyshev  points  and  WENO  ones  use 
20  X  20  grids  of  uniformly  spaced  points.  The  classical  fifth-order  WENO  method,  here  denoted 
WEN05,  uses  the  corresponding  number  of  grid  points  as  if  all  the  subdomains  in  the  Hybrid 
method  were  WENO  subdomains.  The  MR  tolerance  used  is  Cmh  =  5  x  10“^. 

Table  H  shows  that  a  significant  speed  up  is  achieved  when  using  the  Hybrid  method  over 
the  classical  WEN05  for  increasing  resolution.  Eigure  15  shows  the  history  of  the  coverage  of 
the  WENO  subdomains  as  a  percentage  of  the  total  number  of  the  subdomains.  This  percentage 
varies  between  10%  —  20%  for  the  10  X  10  subdomain  partition  and  gets  proportionally  smaller 
by  a  factor  of  2  when  the  number  of  subdomain  partition  increases  by  a  factor  of  2  in  each 
direction. 
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Number  of  subdomains 

Grid  size 

Hybrid  S12W20 

WEN05 

Speedup 

10x10 

200x200 

265 

282 

1.06 

20x20 

400x400 

2009 

2762 

1.37 

40x40 

800x800 

14410 

26090 

1.81 

80x80 

1600x1600 

112900 

253996 

2.24 

Table  II:  CPU  timing  in  seconds  and  speedup  factor  for  the  Shock- Vortex  problem  at  time  t  =  0.6  as 
computed  by  the  Hybrid  (with  constant  Cmr  =  5  x  10“^)  and  the  WEN05  methods. 

Reducing  the  MR  Tolerance  level  of  last  experiment  to  Cmr  =  5  x  10“^,  while  keeping  all  the 
other  parameters  fixed,  leads  to  a  slight  increase  on  the  CPU  time  usage  for  the  Hybrid  method 
as  compared  with  the  previous  case,  as  shown  in  Table  III.  With  a  lower  MR  tolerance,  more 
subdomains  are  classified  as  containing  high  gradients. 


Number  of  subdomains 

Grid  size 

Hybrid  S12W20 

WENO 

Speedup 

10x10 

200x200 

332 

282 

none 

20x20 

400x400 

2239 

2762 

1.23 

40x40 

800x800 

15580 

26090 

1.67 

Table  III:  CPU  timing  in  seconds  and  speedup  factor  for  the  Shock- Vortex  problem  at  time  t  =  0.6  as 
computed  by  the  Hybrid  (constant  Cmr  =  5  x  10“^)  and  WEN05  methods. 


6  Conclusions 

We  have  presented  the  two  dimensional  version  of  the  multi-domain  Hybrid  Spectral- WENO 
method  for  nonlinear  hyperbolic  conservation  laws.  The  rectangular  physical  domain  was  parti¬ 
tioned  into  a  grid  of  subdomains  which  are  either  a  Chebyshev  collocation  grid  for  the  spectral 
method  or  an  uniformly  spaced  grid  for  the  WENO  method.  High  order  Multi-Resolution  Anal¬ 
ysis  is  used  to  measure  the  smoothness  of  the  solution  in  a  given  subdomain  and  a  strategy 
for  switching  the  subdomain  from  spectral  to  WENO  when  the  solution  becomes  non-smooth 
and  vice  versa  was  devised.  The  Hybrid  method  was  tested  with  the  shock-vortex  interaction 
and  the  Richtmyer-Meshkov  instability  (RMI)  two  dimensional  problems.  The  results  were  in 
good  agreement  with  the  classical  fifth  order  WENO  finite  difference  method.  In  these  two  cases 
tested,  the  switching  procedure  devised  performs  well  and  the  coverage  of  the  WENO  method 
were  less  than  20%  of  the  total  number  of  subdomains  resulting  in  a  good  speedup  and  efficient 
method  for  shocked  flow. 

Currently,  we  are  investigating  in  depth  the  role  of  the  Multi-Resolution  Tolerance  Cma  and  its 
effects  on  the  Hybrid  solution.  Furthermore,  the  implementation  of  the  smoothness  measurement 
in  the  spectral  subdomains  can  be  improved  and  is  also  under  investigation.  Since  the  Hybrid 
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method  is  built  in  a  multi-domain  framework,  parallelization  of  the  method  is  fairly  straightfor¬ 
ward;  the  only  required  communications  between  the  subdomains  are  through  a  small  number 
of  ghost  points  used  in  the  WENO  and  Multi-Resolution  Analysis.  A  substantial  speedup  can 
be  achieved  with  a  reasonable  and  balanced  distribution  of  subdomains  to  each  processor. 

We  plan  to  extend  the  Hybrid  method  to  nonlinear  hyperbolic  conservation  laws  systems 
in  three  dimensions  and  with  fifth  and  higher  order  WENO  finite  difference  methods.  One 
application  under  consideration  is  to  study  the  three  dimensional  reactive  flow  over  a  open 
cavity  that  serves  as  the  flameholder  in  a  scramjet  engine.  We  also  plan  to  use  the  Hybrid 
method  for  shocked  flow  over  obstacles  and  other  applications  with  more  complex  geometries. 
Non-conforming  subdomains  with  various  sizes  will  be  investigated  in  this  future  work. 
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(c)  t  =  0.35,  WEN05  Scheme  (d)  t  =  0.60,  WEN05  Scheme 


Figure  9:  Density  p  of  the  Shock- Vortex  interaction  with  Mach  number  Ms  =  3  and  F  =  0.25  at  time  (a) 
t  =  0.35  and  (b)  t  =  0.6  as  computed  by  the  Hybrid  scheme  and  (c)  t  =  0.35  and  (d)  t  =  0.6  as  computed 
by  the  classical  fifth  order  WENO  finite  difference  scheme  (WEN05)  with  1200  x  1200  grid  points. 
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3  4 


(a)  t  =  0.3 


(b)  t  =  0.4 


Figure  10:  Density  p  of  the  Shock- Vortex  interaction  with  Mach  number  Ms  =  6  and  F  =  0.25  at  time 
(a)  t  =  0.3  and  (b)  t  =  0.4  as  computed  by  the  Hybrid  method. 


1)  Rankine-Hugoniot  condition  for  shocks 

2)  Pre-Shock  Temperature  T  =  296  K 

3)  Pre-Shock  Pressure  P  =  0.5  atm 

4)  Xenon  density  pxe  =  2.90  x  10“^ 

5)  Argon  density  pAr  =  0.89  x  10“^ 

6)  Specific  heat  ratio  7  =  | 

7)  Atwood  number  At  =  0.54 

8)  Mach  number  M  =  4.46 

9)  Wave  Length  A  =  3.6  cm 

10)  Amplitude  a  =  1.0  cm 


cnP 


Figure  11:  Initial  Condition  for  the  Richtmyer-Meshkov  Instability  simulation. 
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Regions  of  Interest  : 

1)  Reflected  shock  generated  by  the  shock 

refraction; 

2)  The  penetration  of  the  heavy  (Xe)  to  light 

(Ar)  gas  forms  the  Spike; 

3)  Triple  point  on  the  transmitted  shock; 

4)  A  small  jet  and  Its  vortical  structure; 

5)  The  penetration  of  the  light  (Ar)  to  heavy 

(Xe)  gas  forms  the  Bubble; 

6)  Vortical  rollups  of  the  gaseous  interface. 


Figure  12:  Typical  regions  of  interest  for  the  simulation  of  the  RMI  at  time  t  =  50  /is.  Only  the  lower 
half  of  the  interface  is  shown. 
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(a)  t  =  12.5 


(b)  t  =  25.0  fis 


(c)  t  =  37.5  fis 


(d)  t  =  50.0 


Figure  13:  Contour  plot  of  the  density  p  with  Ms  =  4.46,  A  =  3.6  cm,  a  =  1  cm  at  time  (a)  t  =  12.5  p,s, 
(b)  t  =  25.0  ps,  (c)  t  =  37.5  ps  and  (d)  t  =  50.0  ps  of  the  Richtmyer-Meshkov  Instability  as  computed  by 
the  Hybrid  scheme. 
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(a)  t  =  18.5  iJ,s 


(b)  t  =  52.0  fis 


10  12  14  16  18  20 


(a)  t  =  100.0  jis 


(b)  t  =  100.0  fj,s  with  zoom 


Figure  14:  Contour  plot  of  the  density  p  with  Mg  =  8,X  —  3.6  cm,  a  =  1  cm  at  time  (a)  t  =  18.5  ps,  (b) 
t  =  52.0  ps,  (c)  t  =  100.0  and  (d)  t  =  100.0 /rs  with  zoom  of  the  Richtmyer-Meshkov  Instability  as 
computed  by  the  Hybrid  scheme. 
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Figure  15:  The  time  history  of  the  coverage  of  the  WENO  subdomains  as  a  percentage  of  total  number  of 
subdomains  for  the  four  subdomain  partitions:  (from  left  to  right)  (10  x  10),  (20  x  20),  (40  x  40),  (80  x  80). 


